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ABSTRACT 

We present the first three-dimensional radiation-magnetohydrodynamic simulations of the 
photoionisation of a dense, magnetised molecular globule by an external source of ultraviolet 
radiation. We find that, for the case of a strong ionising field, significant deviations from the 
non-magnetic evolution are seen when the initial magnetic field threading the globule has an 
associated magnetic pressure that is greater than one hundred times the gas pressure. In such 
a strong-field case, the photoevaporating globule will adopt a flattened or "curled up" shape, 
depending on the initial field orientation, and magnetic confinement of the ionised photoevapo- 
ration flow can lead to recombination and subsequent fragmentation during advanced stages of 
the globule evolution. We find suggestive evidence that such magnetic eff'ects may be important 
in the formation of bright, bar-like emission features in H ii regions. We include simple but 
realistic fits to heating and cooling rates in the neutral and molecular gas in the vicinity of a 
high-mass star cluster and show that the frequently used isothermal approximation can lead 
to an overestimate of the importance of gravitational instability in the radiatively imploded 
globule. For globules within 2 parsecs of a high-mass star cluster, we find that heating by stellar 
X rays prevents the molecular gas from cooling below 50 K. 

Key words: H 11 regions - ISM: globules - magnetohydrodynamics - star formation 



1 INTRODUCTION 

When ionising radiation propagates into a non-uniform medium, the 
resultant structure and evolution is very different from the textbook 
development of a spherical Hii region ( |Str6mgren|1939[|Dyson &| 
|Williams|1997[ ). Dense condensations in the neutral gas will slow 
down the propagation of the ionisation front in some directions, 
causing it to "wrap around" the obstacle (see Henney 20()7| § 2.4, 
for a summary). A simple idealisation of this phenomenon is the 
model of a photoevaporating globule: a dense, isolated cloud of 
neutral/molecular gas that is illuminated from one side by a source 
of ionising radiation. The earliest analytical studies of globule pho- 
toevaporation ( |Oort & Spitzer|1955[[Dyson|1968[[Kahn|1969[ l al- 
ready identified the three principal aspects of their evolution: (1) the 
radiation-driven implosion of the neutral/molecular gas, possibly 
giving rise to enhanced star formation; (2) the transonic flow of 
ionised gas back towards the radiation source, capable of driving 
turbulence in the H ii region; (3) the acceleration of the residual neu- 



* Based in part on numerical simulations carried out using the Kan Balam 
supercomputer, operated by the Depailamento de Supercomputo, Direccion 
General de Servicios de Computo Academico, Universidad Nacional 
Autonoma de Mexico. 



tral globule away from the radiation source due to the back reaction 
of the ionised photoevaporation flow (rocket effect). Later semi- 
analytical and numerical wor k ([Sandford et al.|1982[ Tenorio-Tagle| 



|& Bedijn|1981[|Bertoldi|1989||Bertoldi & McKee|1990| l investigated 
the globule evolution in greater detail and exhaustively explored the 
two-dimensional parameter space of initial globule column density 
versus ionisation parameter (ratio of ionising photon density to par- 
ticle density). More recent studies have concentrated on the effects 
of the diffuse radiation field and multiple ionising sources l |Cant6| 
|et al.|1998„Pavlakis et al.|2001[|Cerqueira et al.|2006| ), together with 
detailed studies of the evolution of non-uniform globules, applied to 
gravitational collapse and star formation fKessel-Deynet & Burkert| 
[Gonzalez et al.|2 005 , Esquivel & Raga 2007 ; Motoy ama et al.| 



2003 



2007 



The principal application of globule models has been to the 
study of structures such as bright rims, pillars, and "elephant trunks", 
which are frequently found at the periphery of H n regions around 



young massive stars (Pottasch'1956 


Bedijn & Tenorio-Tagle|1984| 


ICarlqvist et al.|1998; Williams et al. 


2001 1. Two principal modes of 



formation have been considered for the globules: they may simply 
be pre-existing dense cores in the molecular cloud l ,Reipurth_1983j 
[Mellema et al.|2006a^ , or they may be the result of hydrodynamical 
instabilities of the ionisation front and preceding shocked neutral 
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shell (|Spitzer| 1 954| | AxfordI 1 964| |Giuliani| 1 9791 IGarcfa-Segura& 



Francol 1996[|Williams|200 2 , Mizuta et al. 2006 , Whalen & Norman 



2008[ l. The equilibrium globule structure is found to be almost 



identical in the two cases ( Williams et al.|20^T ]l. A second domain 
of application for globule models is to the dense knots seen in some 
planetary nebulae (|Mellema et al.|1998"|[L6pez-Martm et al.|2001[ 



|O'Delletal.|2002||2005t 

Although some previous studies have considered the effects 
of magnetic fields on globule evolution, most such work has been 
qualitative and heuristic only jBertol di 1989, Carlq vist et al.|2003| 
|Ryutov et al.||2005] >. The only published numerical calculations 
have been carried out in one or two dimensions ( [Williams & Dyson] 
|2001 [ [Williams|2007[ l, which is not sufficient to capture the mag- 
netic field geometry in the general case. One-dimensional studies 
of magnetohydrodynamic ( MHD) ionisation fronts jRedman et al.| 
[T998i[Wimams et aL|2000[|Winiams & Dyson|2001t have shown 
that their jump conditions and internal structure can be considerably 
modified from the pure hydrodynamic case, while a recent three- 
dimensional simulation of the evolution of an MHD H n region in a 
uniform medium ( [Krumholz et al.|2007^ has shown that magnetic 
fields dominate the dynamics at late times. 

In this paper, we present the first three-dimensional radiation- 
MHD simulations of the evolution of a photoionised globule. In 
§[2] we describe the numerical algorithms and initial conditions that 
we employ. In §[5] we describe in detail the results of the simula- 
tions with different magnetic field strengths and orientations. In §|4] 
we compare our simulations with observations of photoevaporated 
globules and similar objects. 



equation for hydrogen ionisation/recombination: 



dt 



+ V • (nj) 



,Lc{T)+ f 



rip fie a{T) - ;je 



hydrogen, and electrons, respectively. Additionally, a(T) and C(T) 
are respectively the radiative recombination and collisional ionisa- 
tion coefficients, while cr,, is the photoionization cross-section and 
Jy is the local mean intensity of the ionising radiation field, both 
functions of the photon frequency v. The direct contribution of a 
single, point-like radiation source, of luminosity X*, and located at 
r,, to the local radiation field at a point r is given by 



An\r- r.p 



with 



: 

Jo 



''n('^> + -'^^r) 0"v ds. 



(6) 



(7) 



where e,. is the unit vector (r - r,)l\r - r,\ and s is the distance 
along the straight-line path between f, and The diffuse field due 
to ground-state recombinations is treated in the standard on-the- 
spot approximation ( jOsterbrock & Ferland|2006) , in which it is not 
explicitly included in i,, and the case-B value for a{T) is used. 



2.2 Implementation 



2 NUMERICAL SIMULATIONS 



We carry out our simulations using a new code Phab-C^, which cou- 
ples an Eulerian Godunov MHD code ( |De Colle & Rag a 2005 lOQ^ 
with the radiative transfer/photoionisation code C^-ray ( ,Mellema| 
let al.|2006b|al ). 



2.1 Basic equations 

The equations solved are schematically as follows: 



dp 
dt 



+ V • {pv) = 



dpv 
Ik 



de 
dt 



+ V-((e + p„)fl'-(0'-S)fi) = //-L 



dB 
Ik 



+ W ■{vB-Bv)=0 
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(4) 



where p is the mass density, it is the velocity vector, pta, = p^^^ + B^ 1^ 
is the (magnetic + thermal) total pressure, / is the identity matrix, B 
is the magnetic field (in units of Gauss/y4;r), e is the total energy 



defined as e 



^pv^ + {Bi^ (with 7 = 5/3), and L and 



y_l/'gas ' 2^^" ' 2 

H are respectively the microphysical cooling and heating rates, 
which are functions of the local gas and radiation conditions. The 
above equations represent the conservation of mass (jTJ, momentum 
(|2]l, energy ([3| and magnetic flux (j4j. They are combined with an 



All calculations are performed on a fixed, regular Cartesian grid in 
two or three dimensions. The MHD Riemann solver uses a standard 
second-order Runge-Kutta method for the time integration and a 
spatially second-order reconstruction of the primitive variables at 
the interfaces (except in shocks). Lapidus viscosity is applied to 
fluxes at cell interfaces, following [Colella & Woodward) p"984[ l. 
The constrained transport (CT) method (e.g., |T6th|2000^ is used to 
conserve V • 6 = to machine accuracy. 

The C^-ray (Conservative-Causal ray, Mellema et al.|2006b 
algorithm uses a short-characteristic method to calculate the ra- 
diative transfer of ionising radiation (equation (^), together with 
an explicitly photon-conserving iterative technique for solving the 
ionisation rate equation (equation (|5]l), which allows one to use 
timesteps that are much longer than the ionisation/recombination 
timescales without loss of accuracy ( ,Iliev et al.|2006^ . 

The radiative heating and cooling terms {H and L in equa- 
tion ([3j) are calculated explicitly, using a local fractional timestep 
technique to accurately treat those cells undergoing rapid tempera- 
ture change. Heating is principally due to absorption of stellar radi- 
ation by gas or dust: ionizing extreme ultraviolet (EUV) radiation 
in the ionized region, non-ionizing far ultraviolet (FUV) radiation 
in the neutral gas and x rays in the most shielded molecular gas. 
Cooling is principally due to collisionally excited line radiation of 
ions, atoms, or molecules. All of the various processes that con- 
tribute to H and L, together with their implementation in the code, 
are discussed in detail in Appendix|A] The individual components 
of the Phab-C^ code have been extensively tested against standard 
problems (De Colle,2005||De Colle & Raga|2004, Iliev et al. 2006} . 
In addition, Arthur et al.| ( 2009 1 present further test cases for the 
combined algorithm. 
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Figure 1. Initial conditions for the globule simulations. The grayscale shows 
the initial gas density in the midplane z = 0.5 pc. 



2.3 Initial and boundary conditions 

The initial setup for our calculations is illustrated in Figure[T] Our 
simulations are calculated in a box of dimensions (x,nax, ymax, Zmax) = 
(2, 1, 1) pc. The initial conditions consist of a spherical neutral 
globule with peak density no = W* cmT^ and characteristic radius 
Tg = 0.2 pc, located at a position rg = (xg, yg, zo) = (0.5, 0.5, 0.5) pc. 
An ultraviolet radiation source with ionising photon luminosity 
Qh = 10*' s"' and effective temperature T^ff = 40, 000 K is located 
at a position r, = (x,,y,,z,) = (0.0,0.5,0.5) pc, so that the initial 
distance of the globule centre from the source is Do = 0.5 pc. The 
stellar FUV photon luminosity is 2fuv = 6.25 x lO'"* s"' and the 
stellar x-ray luminosity is assumed to be Lxj = 2.29 x 10^^'' erg s"' 
with temperature Tx,i = 1.6 x 10' K. In addition, we include a 
distributed x-ray source, corresponding to young low-mass stars in 
an accompanying star cluster, with a luminosity of Lx^ = 3.47 x 
10^^^^ erg s"' and Tx,2 = 2.5 x 10' K, with the x-ray flux from this 
cluster calculated assuming a softening length of 0.3 pc. These x- 
ray properties were chosen to be consistent with observations of 
the Trapezium cluster in the Orion Nebula ( [Flaccomio et al.|2003[ 
[Feigelson et al.|2005[[Ste"lzer et al.|2005t . 

The globule has a smooth Gaussian density profile and is sur- 
rounded by a uniform ambient medium of density = 100 cm"^, 
such that the initial density as a function of position r in the box is 
given by 

n(r) = n,, + (no - "J exp(-|r - fol'/rl). (8) 

The initial ionisation fraction is set to the constant value of 0.01, 
while the initial temperature is set inversely proportional to the 
density, so as to give a constant pressure, with a temperature in the 
ambient medium of T.^ = 1500 K and a minimum temperature in the 
globule of To = 15 K. The mean mass per hydrogen nucleon is set 
torn = 1.3mp = 2.174 X lO^^"* g. 

In this initial investigation of the problem, we vary only 
those parameters describing the magnetic field, keeping all other 
globule properties fixed. The globule mass for the adopted Gaus- 
sian density profile is Mg = n^^-r^mng = 14.3 Mq, whereas the 
mass of the diffuse ambient medium in our simulation box is 

Ma = XmaxS^maxZmaxna"! - n.^M^jllo = 6.3 Mq. 

We assume that the magnetic field is initially uniform, with 
magnitude Bq, lying in the xy plane at an angle Q(, to the x-axis 
(direction of ionising photons). We consider three strengths for the 
magnetic field: zero field (So = 0), weak field (Bq = 59 yuG), and 
strong field (Bo =186 ^G). The strong field is similar to the field 
strengths that have been observationally inferred in dark globules 
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Table 1. Parameters for simulation runs 



Run 


Bn (uG) 




fo(°) 


Dim 




S80H 


186 


0.01 


80 


3 


510 


S80L 


186 


0.01 


80 


3 


254 


SOOL 


186 


0.01 





3 


254 


S45L 


186 


0.01 


45 


3 


254 


W80L 


59 


0.1 


80 


3 


254 


ZOOL 





DO 





3 


254 


S80S 


186 


0.01 


80 


2 


1002 




Figure 2. Idealised sketch of a photoevaporating magnetic globule. Ionised 
gas is shown in red, neutral/molecular gas is shown in green, magnetic field 
lines are purple. The sketch shows the configuration for an initial magnetic 
field orientation of 9o ~ 80° , but qualitatively similar structures result from 
other field orientations. 

(e.g., [Wolf et al.|2003| l. The plasma /? parameter (ratio of gas pressure 
to magnetic pressure) is initially constant throughout the simulation 
box, being y6o = 0-1 in the weak field case and /?o = 0.01 in the 
strong field case. Our most detailed investigations are of the almost 
perpendicular field case (So = 80°), but we also present simulations 
with (?o = 0° and So = 45°. The three-dimensional simulations are 
carried out at resolutions of up to 510 x 255 x 255 cells, and we 
also carry out one simulation in two-dimensional slab geometry 
at a resolution of 1002 x 501 cells. The parameters for each run 
are summarised in Table [T] We do not present simulations with 
(So = 90°) since the exact symmetry about the y = Q plane causes 
numerical difficulties at our lowest spatial resolution due to the 
formation of a thin dense sheet that is exactly aligned with the grid 
axes. 

For all of these runs, standard outflow boundary conditions are 
adopted on all faces of the simulation cube. The effects of varying 
these boundary conditions are discussed in § |3.6| 



3 RESULTS 

In the following description of the simulation results, particularly 
in figure captions, we will use a vocabulary of directions relative to 
the position of the globule, which should be interpreted as follows. 
"Above" and "below" refer to larger and smaller values of y, respec- 
tively. "Behind" and "in front" refer to larger and smaller values of 
X, respectively, so the ionising star is "in front" of the globule. "To 
the side" means away from the z = 0.5 pc symmetry plane. 

Colour images of the optical appearance of the simulations are 
calculated as in |Mellema et al.| l f2006a^ , with each of the three colour 
channels representing the surface brightness in a different emission 
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line: [Nii] 6584 A (red). Ho- 6563 A (green), and [Oiii] 5007 A 
(blue). Since our simulations do not explicitly follow the ionisation 
state of elements other than hydrogen, the ion fractions of N and 
O were approximated as fixed functions of the hydrogen ionisation 
fraction (Henney et al. 2005). The effects of dust absorption are 
included in calculating the images, asuming Orion Nebula dust 
properties ( [Baldwin et al.|199T) . A sketch of a typical stage from 
our simulations is shown in Figure |2] which explains some of the 
terms used in the description of our results. 



3.1 Zero field: ZOOL 

For reference, we first briefly discuss the globule evolution in the 
absence of a magnetic field (run ZOOL), in which case, the globule 
evolution is cylindrically symmetric. Various snapshots of the evo- 
lution are shown in Figure |3] Initially, the low density halo of the 
globule is ionised by a fast-moving R-type ionisation front, which 
transitions to D-type after about 1000 years, driving a convergent 
shock into the globule. The shock passes through the globule after 
about 50,000 years, after which time of maximum compression 
(second panel in Fig. [3j the globule is coherently accelerated by 
the back reaction of the transonic photoevaporated flow that leaves 
the ionisation front (rocket effect). At the same time, the ionised 
photoevaporated flow sweeps all the ambient material off the grid in 
the direction facing the star After the globule shock bounces off the 
symmetry axis, the neutral globule reexpands somewhat as it recedes 
from the star. All these features of the non-magnetic evolution have 



been well-studied by previous workers i Bertoldi 1989 |Lefloch & 



|Lazareff' T994l|Mellema et al.|1998[|Pav akis et a .120011 IWilliams 
[et al.,200ll . 



3.2 Strong perpendicular field: S80H, S80L, S80S 

In the ideal MHD approximation, matter is perfectly coupled to the 
magnetic field lines. Therefore, when the magnetic field in the neu- 
tral globule is very strong, magnetic pressure and tension strongly 
resist any movement of the gas in directions perpendicular to the 
original field direction, whereas no such magnetic support exists 
to oppose gas motions along the field lines. As a result, the implo- 
sion of the globule by the ionisation-driven shock front is highly 
anisotropic: the globule is swiftly flattened along the (/-direction, 
whereas the radius of curvature in the xz plane remains large, giv- 
ing a disk-like structure to the globule. This can be appreciated in 
Figure|4] which shows the structure of the ionised emission from 
run S80H. At the earlier time shown (20,000 years, top row) the 
ionisation front already shows strong departures from cylindrical 
symmetry, such as the "whiskers" seen in the xz plane, while after 
the point of maximum compression (50,000 years, middle row), the 
flattening is considerable (aspect ratio of 5 : 1). 

A further difference from the non-magnetic case is seen in the 
behaviour of the ionised photoevaporation flow. Unlike in run ZOOL, 
the photoevaporation flow cannot sweep all the ambient gas off the 
grid because the ambient magnetic pressure is too high. Instead, it 
forms a shocked shell, which is anchored by the large-scale magnetic 
field that threads the simulation box, as can be best seen as strong 
blue emission at the left-hand side of each panel in the lower rows 
of Figure |4] 

In order to investigate the structure of the globule in greater 
detail, we have carried out a two-dimensional simulation in slab xy 
geometry (run S80S), which allows us to achieve a ~ 2 times finer 
spatial resolution than in the highest resolution three-dimensional 



( 



ZOOL, t = 30,000 years, 6 = 350°, (f) = 350° 



ZOOL, t = 60,000 years, 6» = 350°, cf) = 350° 




ZOOL, t = 90,000 years, = 350°, cf) = 350° 




ZOOL, t = 120,000 years, 6 = 350°, </> = 350° 




ZOOL, t = 150,000 years, = 350°, (/> = 350° 



Figure 3. Sideways view close to the z-axis of the ionised emission from 
run ZOOL for a sequence of evolutionary times from 30,000 to 150,000 years 
(top to bottom). The colours represent the surface brightness in three optical 
emission lines, as discussed in the text: [Nil] 6584 A (red), Htt 6563 A 
(green), and [Oui] 5007 A (blue). Thus, yellow-orange colours trace the 
ionisation front on the surface of the globule, while blue-green colours trace 
highly ionised gas. 
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Figure 4. As Fig.[3]but showing multiple views of ionised emission from the early evolution of run S80H (strong perpendicular field) for 20,000 years (top), 
50,000 years (middle), and 80,000 years (bottom). Left column shows the view from the side, slightly in front and slightly above. Central column shows the 
view from behind and below. Right column shows the view from below. See §[2]for the definition of these in terms of the coordinate axes. 



jWimams et al.|2000[[MTiams|2007> since the flow passes through 
both the Alfven and sound speeds. 

The slow magnetosonic speed is zero in the direction exactly 
perpendicular to the magnetic field lines, so the slow-mode shock 
cannot propagate in the symmetry plane (y = 0.5 pc). Instead, a 
dense ridge is formed from the pincer-like convergence of the neutral 
shell from above and below. The ridge accretes both material and 
horizontally oriented field through fast-mode shocks at its upper and 
lower boundaries. The triply- shocked gas in the ridge has densities 
as high as 2 X 10* cm"', 200 times denser than the initial globule. 
The field lines from above and below the midplane have opposite 
directions, so a current sheet forms in the mid-plane, in which 
magnetic reconnection can occur. 

3.2.2 Accelerated globule phase 

By a time of 60,000 years (third row of Fig.|5j, the initial implosion 
shocks have passed through the entire globule, forming a thin, dense 
sheet of material in the midplane. The neutral globule has a much 
more complicated internal structure than in the pure-hydrodynamic 
case due to the mutual interactions of the multiple shocks that have 
reflected from the midplane. However, the head of the globule has a 
configuration that is qualitatively similar to that seen at the earlier 
time: the magnetic field lines are approximately perpendicular to 
the ionisation front everywhere except for the nose that faces the 
star. The neutral dense ridge that was forming at 40,000 years has 
by now been completely ionised, so that the nose now consists of 
material from the neutral shell that formed on the flanks of the 
original globule. The magnetic field lines in these parts of the shell 
were not bent far from their initial (/-orientationQso the mid-plane 

' The effects of the oblique fast-mode and slow-mode shocks are in opposite 
directions and approximately cancel. 



runs. The use of slab geometry is equivalent to imposing djdz = 
for all quantities. This is a reasonably good approximation to the 
mid-plane (z = 0.5) of the three-dimensional runs with a strong per- 
pendicular field, since the flattened form of the globule means that 
d/dz d/dy. The results are shown in Figure[5]for three different 
times. Except where specifically noted, all the following description 
applies equally to the two-dimensional and three-dimensional run. 



3.2.1 Initial radiation-driven implosion 

At the earliest time of 10,000 years (top row of Fig.|5](, a fast-mode 
shock has started to run ahead of the ionisation front, compressing 
both the gas and magnetic field in the globule. This is followed 
by a slow-mode shock, which at this earliest time can be seen 
just starting to detach from the ionisation front at the sides of the 
globule. The ionised portion of the globule is starting to flow back 
towards the star due to the pressure imbalance (the equilibrium 
temperature in the ionised gas is almost independent of density), but 
the photoevaporation flow is not yet fully developed. 

At the intermediate time of 40,000 years (second row of Fig.|5]l, 
both the fast-mode and slow-mode shock are clearly visible, which 
compress the gas by a factors of ~ 2 and ~ 10, respectively, pro- 
ducing a neutral shell with density ~ 20 times greater than the 
initial peak globule density. In addition, the slow-mode shock is of 
the "switch-off" type, in which the post-shock magnetic field has 
no component in the plane of the shock. The dense neutral shell 
is therefore effectively demagnetised and gas-pressure dominated 
(fi ~ 5). Since the neutral shell is thin, the ionisation front and the 
slow-mode shock are approximately parallel, meaning that the B- 
field is perpendicular to the ionisation front over much of the surface 
of the globule. Such fronts, in which the jump conditions reduce 
to the pure hydrodynamic case, have been dubbed "extra strong" 
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Figure 5. Physical structure of run S80S (strong perpendicular field, two-dimensional slab geometry) at evolutionary times of (top to bottom) 10,000, 40,000, 
60,000 and 100,000 years. Left panels show magnetic field lines (grey lines) and gas velocity (black arrows). Right panels show logarithm of total hydrogen 
number density, n (grey scale, see scale bar) and ionisation fraction (dashed contour at a value of 0.5). 



current sheet has disappeared. Thin-shell instabilities start to grow 
in the shocked sheet. Although these instabilities occur in all our 
simulations, their onset is slightly earlier and their growth is more 
vigorous for the higher resolution runs. 

The ionised photoevaporation flow from the head of the globule 
at this time consists of two zones with distinct properties. The 
flow from the flat nose of the globule is very weakly magnetised 
(P = 10-100), whereas the surrounding sheath flow from the globule 
shoulders is closer to equipartition {fi = 1-3). For comparison, after 
ionisation by an R-type front, the quiescent ambient medium would 
have - 0.2. This structure does not remain constant in time, but 
evolves as the ionisation front progresses through different regions 
of the complex magnetic geometry in the globule head, with each 
zone of the photoevaporation flow slowly peeling back towards the 
sides. In both zones the flow is mildly supersonic and accelerating, 
reaching maximum velocities ^ 30 km s"'. Fine-scale structure in 
the magnetic field at the ionisation front produces "streamers" in 
the photoevaporation flow, which typically have lower-than-average 
field strength, together with higher-than-average density and velocity. 
The whole photoevaporation flow is criss-crossed with weak oblique 
shocks. 



The shadow tail behind the globule is accreting gas along the 
field lines from the ambient medium to either side, due to the higher 
thermal pressure in the ionised gas. This accretion flow, with velocity 
10 km s"', recombines in a broad front as it enters the region 
shielded from direct ionising radiation, and then shocks against the 
thin dense (n ^ 10'*-10' cm"') core of the tail. It should be noted 
that our treatment of the diffuse ionising field is very approximate 
and that this will predominantly effect the structure of the tail region 
( [Canto et al.|1998] >. However, globule simulations that include an 
exact treatment of the diffuse field ( [Raga et al.|2009[ > show that the 
differences are not large in most cases. 

At later times, the instabilities in the shocked globule continue 
to grow. In the globule head, these cause fragments of dense neutral 
gas to separate from the main globule, forming isolated smaller 
globules, as can be seen at at time of 100,000 years in Figure |5] 
In the globule tail, the instabilities cause portions of the shocked 
sheet to exit the zone behind the head that is shadowed from direct 
ionising radiation. When this happens, the tail portion begins to 
be photoevaporated, but this is a transient phenomenon since the 
rocket effect then drives the neutral sheet back into the shadow zone. 
The process repeats when the sheet leaves the other side of the 
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shadow, producing quasi-periodic oscillations on a timescale of tens 
of thousands of years. 

3.2.3 Late-time recombination of the shocked shell 

At times later than 60,000 years, we begin to see significant de- 
viations between the behaviour of the two-dimensional and three- 
dimensional runs0This is due to the evolution of the shocked shell 
that forms by the interaction of the ram-pressure dominated photoe- 
vaporation flow with the magnetically dominated ambient medium. 
In the three-dimensional run, the cavity that is swept out by the 
photoevaporation flow has a finite size in the z direction of about 
0.7 pc, roughly equal to the globule radius of curvature in the xz 
plane, whereas in the two-dimensional run, the cavity is implicitly 
unlimited in its z extent. As the globule is pushed away from the 
star by the rocket effect, the flow cavity moves with it, vacating the 
region with x < 0.5 pc. In the three-dimensional run, the ambient 
medium that was pushed aside by the photoevaporation flow then 
closes in behind the cavity like curtains, driven by the magnetic pres- 
sure and tension in the cavity walls. At a time of ^ 80, 000 years, 
the two sides of the shocked shell collide with one another in the 
symmetry plane z = 0.5 pc at a relative velocity of ^ 40 km s"', 
forming a dense ribbon (n ^ 1500 cm"') of ionised gas, elongated 
along the y-axis. This can be seen in the bottom row of Figureplfor 
the high-resolution run S80II and in the top panel of Figure |6[f or 
the low-resolution run S80L. 

A few thousand years later, the density in the shocked ribbon 
is sufficiently high to trap the ionisation front, causing parts of the 
ribbon to recombine, as do sections of the cavity wall. After this 
time, the central part of the original globule is in the ionisation 
shadow of the ribbon, so this too recombines, the photoevaporation 
flow ceases, and the now-overpressured neutral globule begins to 
expand and dissipate. This can be seen in the second row of Figure|6] 
However, by a time of 120,000 years (third row of Fig.|6]l, all of this 
recombined material has been reionised and photo-evaporated, so 
that the original globule is once again illuminated by the ionising 
radiation. At the same time, as with the 2D models discussed above, 
the globule's trajectory begins to veer increasingly away from the 
X-axis as its motion becomes closer to parallel to the magnetic field 
lines. The complex flow structure in this late stage is seen in the 
bottom row of Figure [6] 

3.2.4 Thermal behavior of the globule 

Figure|7]shows the evolution with time of the gas density and tem- 
perature for model S80L. The initial conditions are shown in the 
top-centre panel, appearing as a diagonal line corresponding to 
nT = constant, since the initial configuration is at constant pressure. 
This initial configuration is not in thermal equilibrium, mostly lying 
above the region delimited by the PDR equilibrium curves (thin 
solid lines in the figure), and with a cooling time (dotted lines in the 
figure) of a few thousand years. 

By a time of 10,000 years (top-right panel), the initial config- 
uration has been modified in various ways. A portion of the gas 
has been ionised, giving the horizontal band at roughly constant 
temperature (logT ^ 3.9). Another portion of the gas lies in the 
thin cooling zone behind the convergent fast-mode shock that is 
driven into the globule, giving the curved band of material between 

^ We have only been able to fully investigate this late-time evolution in the 
lower resolution run S80L. 



(log n, log r) = (3.0,2.0) and (4.1,2.6). The small amount of ma- 
terial around (log n, log T) = (4.0,3.1) has been shocked by the 
slow-mode shock that at this time has just detached from the ionisa- 
tion front on the flanks of the globule (see Fig.[5j. The remainder 
of the globule gas (logn = 2.5 to 4.0) has reached its equilibrium 
PDR temperature of log T = 1.5 to 2.3, depending principally on 
the extinguishing column, which at this early stage of the implosion 
is still moderate (Ay - 3 to the center of the globule). 

By a time of 50,000 years (bottom-left panel), the core of the 
globule reaches its maximum compression, having been reshocked 
to densities ~ 10* cm"'. Due to the 4 times lower linear resolution, 
this is about a factor of 3 lower than the maximum density seen 
in the slab-symmetric two-dimensional model S80S discussed in 
§ |3.2.1| The density-temperature distribution of model S80S at 
t = 50, 000 years is also shown in Figure |7] (bottom-centre panel). 
It is very similar to S80L, except for extending to slightly higher 
densities in the imploded globule. The greater noise apparent in the 
S80S graph is due to the fact that, despite the higher resolution, it 
has considerably fewer grid cells than S80L. 

Our simulations do not include self-gravity, but an estimate 
of whether gravitational efi'ects may be important can be ob- 
tained by considering the leans instability criterion and the free- 
fall time. The heavy solid lines in Figure [7] show the threshold 
for leans instability for globules of 10 Mq and 1 M^. Mj = 
2Mo(c70.2 kms"')'(;j/10' cm"')""^. This can only be taken as 
indicative, since it does not include the elfects of magnetic fields nor 
turbulent motions, which would help suppress the instability. The 
total neutral mass of the globule at this point is about 6 Mq, but only 
a fraction of this mass has the highest densities. Furthermore, the 
dense gas is not in a compact configuration, but rather is in the form 
of a ridge, which would further reduce its effective mean density for 
the purposes of triggering the instability, so it seems likely that the 
globule would resist global gravitational collapse. An additional ar- 
gument is that the free-fall time (oc n"''-) is of order 100, 000 years, 
which is long compared with the dynamical evolution time. It may 
be that a fraction of the globule mass is compressed to such high den- 
sities that it undergoes prompt gravitational collapse, but addressing 
this question would require a much higher spatial resolution than 
we can achieve in our fixed-grid simulations. 

By a time of 110,000 years (bottom-right panel), the glob- 
ule has rebounded from its maximum compression and entered its 
equilibrium acceleration phase (§ |3.2.2[ l. At this point, the density- 
temperature distribution consists of a horizontal band of ionised gas 
at logT ^ 3.9, from which descend 3 strands of neutral material 
towards lower temperatures and higher densities. The rightmost 
strand corresponds to the bulk of the original globule material, the 
central strand corresponds to the shocked tail, whereas the leftmost 
strand, which is very faint since it contains little mass, corresponds 
to the partially ionised transonic accretion flow onto the tail. 

It is obvious from Figure |7] that the widely used assumption 
of an isothermal equation of state for the neutral/molecular gas is a 
very poor approximation for much of the evolution of the globule. 
The combined effect of heating by shocks and stellar x rays keeps 
the bulk of the dense gas (> 10' cm"') at temperatures in the range 
50 to 100 K, whereas less dense neutral gas (~ 10"* cm"') is fre- 
quently heated to 100 to 500 K by FUV radiation. This will have 
important consequences for questions of gravitational instability and 
triggered star formation. For instance, Esquivel & Raga (2007) as- 
sumed a constant T = 10 K for the neutral gas in simulations of the 
photoevaporation of self-gravitating clumps, finding that the Jeans 
instability was triggered in some cases. If they had used T = 50 K, 
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Figure 6. As Fig.|4]but showing multiple views (left, centre and right columns) of the late-time evolution of a run S80L (strong, nearly perpendicular field) at 
times from 65,000 (top row) to 140,000 years (bottom row). The first row is for a time shortly after the last time shown in Fig.[4]and illustrates a weaker growth 
of instabilities in the shocked globule in this lower resolution run. The second row shows the moment in which material in the dense shocked ribbon and cavity 
walls have recombined to shadow sections of the original globule (see text), whereas, by the time of the third row, this material has been reionised. At the time of 
the fourth row, the globule is about to leave the top of the grid after being deflected by the magnetic field. 



then the Jeans mass would have been 5^'" = 1 1.2 times higher and 
their clumps would have remained gravitationally stable. 

One shortcoming of our simulations is that we use one- 
dimensional models to calibrate the heating and cooling functions 
(Appendix [a), whereas in reality multidimensional effects may be 
important. In particular, during its initial compression phase the 
neutral globule is compressed to a thin sheet with a high optical 
depth to stellar photons, which propagate parallel to the sheet, but a 
much lower optical depth in the perpendicular direction. This may 
enhance the cooling efficiency in the dense gas by allowing optical 
and near-infrared radiation to escape more easily. However, in prac- 
tice this effect is likely to be small since the cooling for T < 1000 K 
is dominated by far-infrared and millimeter lines, for which the 
sheet is optically thin in all directions. 



3.3 Effects of varying the initial field angle 

In order to investigate the effect of the initial magnetic field geometry 
on the globule evolution, we have carried out runs with 0q =45° 
and 0°. These are both at the lower resolution of 255 x 127 x 127. 



3.3.1 Strong parallel field: SOOL 

When 6q = 0°, the x-axis, (j/, z) = (0.5, 0.5) pc, becomes a symmetry 
axis of the problem, as in the non-magnetic case, so that the globule 
remains cylindrically symmetric throughout its evolution. One im- 
portant difference with respect to the non-magnetic case is that the 
support provided by the magnetic field opposes the lateral compres- 
sion of the neutral globule, producing a broader, snubber globule 
head. The internal structure of the compressed globule is much sim- 
pler than in the perpendicular field case: consisting of a thin shell 
behind the ionisation front, with peak density « 10* cm"^, plus 
a lower density core with n ^ 5000 cm"'. The shell is moderately 
magnetically dominated (fi ^ 0.3), while the core and the tail (which 
maintains its ambient density of 100 cm"') are overwhelmingly mag- 
netically dominated (/3 ^ 0.001). The field lines within the globule 
never get bent through large angles and there is no equivalent of the 
current sheet that forms in the Bg = 80° runs. There is no accretion 
flow into the tail because the thermal overpressure of the ionised 
gas at the sides of the globule is not sufficient to laterally compress 
the tail's longitudinal magnetic field. The ionised emission from the 
SOOL run for various evolutionary times is shown in the left column 
of Figure [8] 

The magnetic field in the shocked neutral shell is approximately 
perpendicular to the globule surface over the entire globule head, so 
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Figure 7. Evolution of thermodynamic quantities. The grayscale images show the fraction of the total gas mass in the simulation that lies within a given 
logarithmic interval of density (.v axis) and temperature (y axis). Results are shown for a sequence of evolutionary times of model S80L. Thin solid lines show 
equilibrium temperatures for neutral/molecular gas, at different distances, r in parsecs, and extinguishing columns. Ay in magnitudes, from the ionising star — top 
to bottom: {Av,r) = (0,0.3), (0, 1.0), (2,0.3), (2, 1.0), (10,0.3), and (10, 1.0). Dotted lines show contours of radiative cooling time — top-right to bottom-left: 
/cool = 100, 1000, 10, 000, and 100, 000 years. Heavy diagonal solid lines indicate Jeans masses of 10 Mq (upper line) and 1 M© (lower line) — Jeans instability 
requires a clump of the requisite mass to be below the line. The region to the right of the vertical dashed line corresponds to a gravitational free-fall timescale 
< 100, 000 years. 



that ionisation front is again of the "extra-strong" type. The ambient 
magnetic field channels the photoevaporation flow towards the sym- 
metry axis by means of a focusing shock. At around 60,000 years 
the sides of the shocked shell become optically thick to ionising radi- 
ation (shown in the second row of Fig.[8]l, and a few thousand years 
later so does a dense ionised knot that has formed on the symmetry 
axis. This causes the recombination of the photoevaporation flow 
(third row of Fig.[8]l in a similar way to the 6o = 80° runs (§ |3.2.3^ , 
although the recombination period lasts much longer in this case 
(bottom two rows of Fig.[8j. 

3.3.2 Strong slanted field: S45L 

The evolution of the 0q = 45° run shares elements of both the 
9q = 80° (§ |3.2| > and 6o = 0° runs (§ |3.3.1| l, but also has unique 
features of its own. The ionised emission from the S45L run is 
shown in the right column of Figure[8]for a sequence of evolutionary 
times. The problem no longer has an xz symmetry plane and the 
globule acceleration is not strictly radial away from the ionising 
star. The magnetic field is not initially strong enough to force the 
globule to move like a bead on a wire, but magnetic tension forces 



are sufficient to deflect the globule motion by =^ 15° from the x-axis. 
A feature analogous to the dense ridge that formed in the 9q = 80° 
runs at the symmetry plane now forms at the top corner of the 
globule head, where the initial magnetic field is tangential to the 
globule surface. However, unlike in the 9o = 80° case, the ridge does 
not have a large radius of curvature in the xz plane and the globule 
does not become significantly flattened. Instead, the unbalanced 
forces due to the one-sided ionising illumination of the ridge cause 
it to "curl up" away from the ionising star, so that the magnetic 
geometry of the head becomes similar to that in the do = 0° run, 
with the photoevaporation flow streaming off the ionisation front 
along approximately radially oriented field lines. As with the other 
strong field runs, the shocked ionised shell eventually traps the 
ionisation front, leading to recombination of the photoevaporation 
flow and a complex subsequent evolution. 



3.4 Weak perpendicular field: W80L 

This run has an initial magnetic field that is -^10 times smaller 
than in the strong-field runs, giving an initial plasma parameter of 
po = 0.1. Although the magnetic pressure dominates the thermal 
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SOOL, t = 120,000 years, fl = 350°, ^ = 350° S45L, t = 120,000 years. 6 = 350°, = 350" 



Figure 8. As Fig.[3]but showing views of ionised emission from runs SOOL 
(strong parallel field, left column) and S45L (strong slanted field, right 
column), for a sequence of evolutionary times. View directions are the same 
as in the leftmost column of Fig.|4] 

pressure in the initial globule, the thermal pressure of the ionised 
photoevaporation flow is much greater than this and so the magnetic 
field plays a minor role in the subsequent evolution of the globule, 
which is more similar to the non-magnetic case (§ |3.1[ l than to the 
strong-field case (§ |3.2[ l. In particular, the flattening of the globule 
along the field lines (§ |3.2.1[ l is much less extreme with an aspect 
ratio of about 2:1. In addition, no current sheet forms at the sym- 
metry plane (§ |3.2.1[ l, and the photoevaporation flow is powerful 
enough to drive all ambient material off the grid, so there is no 
recombining shocked shell (§ |3.2.3^ . The evolution of this case is 
shown in Figure |9] 

3.5 Comparative study of global globule evolution 

Figures [To] to [13] compare the evolution of various global average 
properties between the different runs. For all the runs, the mass of 
neutral gas (Fig. [TO] top) falls in an identical manner during the 
initial stages of shock compression, as the globule begins to accel- 
erate (Fig.|l 1[>, but after about 20,000 yr small differences become 
apparent. The strongly magnetised globules all evaporate faster than 
the weakly magnetised and unmagnetised globules, but for different 
reasons. In the 9o = 0° case (run SOOL), where the effect is the 
strongest, it is because the lateral compression of the globule is 
retarded by the magnetic field (§ |3.3.1^ , so that the globule presents 
a larger cross-section to the ionising radiation, increasing the pho- 
toevaporation rate. In the Oq = 80° case (runs S80L and S80H), the 
evaporation rate does not begin to exceed the non-magnetic case 
until after 40,000 years, at which point the globule center-of-mass 
velocity has stopped increasing (Fig.[TT] top) due to braking by the 
threaded field lines, so that at a given time the globule is closer to 
the star than in the non-magnetic case and is therefore exposed to a 
greater ionising flux. 

An extended plateau is seen in the mass of ionised gas (Fig.|10[ 
bottom) after ^ 30, 000 years, during which time the globule is in 



its equilibrium cometary phase. The amount of ionised gas that is 
retained on the grid is much higher in the strongly magnetised runs, 
due to the trapping of the ionised photoevaporation flow in a shell 
(§ |3.2[ l. Apart from in the Og = 80° strong field case, where magnetic 
braking is efficient, all the runs show a very similar evolution in the 
globule center-of-mass velocity along the x-axis (Fig.[TT] top). In the 
slanted field runs (S45L, S80L), the globule develops a lateral veloc- 
ity along the y-axis (Fig.^^ bottom). Interestingly, the magnitude 
of the lateral velocity is similar in the two cases, although the time 
evolution is somewhat different. For 9o = 45° the globule moves 
approximately in a straight line with constant (Vj)/(Vv), whereas 
for Oq = 80°, (V„)/(V.r) increases with time as the globule's path 
curves away from the radial direction. 

After a time between 60 and 80,000 years, depending on the 
field inclination Oq, the equilibrium cometary phase for the strongly 
magnetised globules is modified due to recombination in the shocked 
ionised shell. In the cases where this effect is greatest (S45L and 
SOOL), this can be seen as an increase in the neutral mass (Fig.|10[ 
top), which in the case of do = 0° shows multiple epsisodes of 
recombination and reionisation. 

The mean magnetic field for the ionised or neutral component 
is calculated as (|6,|) = j {Bjlxj dV/ j xj dV, where Bj is 5^. or 
and Xj is the ionised or neutral fraction of hydrogen, as appropriate. 
The evolution of the mean magnetic field in the neutral gas is shown 
in Figure [T2] where it can be seen that in the strong field cases 
the variations from the initial value (Table[T]l are rather slight (10- 
20%)|^ During the initial implosion phase, there is competition 
between the leading fast-mode shock, which amplifies 6, and the 
following slow-mode shock, which attenuates B. The result is a small 
gradual increase in the mean B. The first rebound after maximum 
compression reverses this trend, but the late-time evolution in the 
strong-field runs is complicated by the recombination of the ionised 
shell. In the weak field case (W80L), the field amplification is more 
significant, reaching a factor of two. 

Figure^jshows the evolution of the mean plasma fl parameter, 
(/?) = <-Pgas)/(^mag) for both the neutral gas (top) and ionised gas 
(bottom), normalised by the initial value (JUg = 0.1 for the weak- 
field runs, j8o = 0.01 for the strong-field runs). In the ionised gas, 
{pyiPo is initially in the range 40-60, falling to 10-20 at later stages. 
Thus, on average, the magnetic pressure exceeds the gas pressure in 
the ionised region by a factor of 2-10. However, this represents a 
compromise between the photoevaporation flow from the globule 
head, where the gas pressure dominate^ (^0 > 1), and the flanks 
of the tail, where the magnetic field dominates {f} < 1). In the 
weak field case, the ionised (Ji)/fio falls off somewhat during the 
equilibrium cometary phase from the initial value of =^ 40, but the 
magnetic contribution to the total pressure in the ionised gas remains 
small. 

In the neutral gas (Fig.^]top), (fi)//io tends to increase during 
the globule evolution, but in the strong field cases the magnetic 
pressure remains dominant over the thermal pressure by about one 
order of magnitude. During the initial implosion, the increase in 
(y6) is due to pressurisation of the globule by shocks. The strong, 
almost perpendicular field run (S80L) shows a larger boost in the 



' In part, this is due to the fact that much of the contribution to the mean 
field comes from the globule tail, where the magnetic field is relatively 
unperturbed. 

* Since the flow is supersonic, the gas pressure is in turn dominated by the 
ram pressure. 
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Figure 9. As Fig.|3] but showing multiple views of run W80L for a sequence of evolutionary times between 30,000 and 120,000 years. View directions are as in 
Fig.|4] 



neutral (fi) than the other models, which is probably due in part to 
reconnection in the current sheet that forms in the globule midplane. 

Although the evolution of many of the global properties in the 
strong perpendicular field case is almost indistinguishable between 
the low-resolution (S80L) and high-resolution (S80H) runs, there 
are small but significant differences after the point of maximum 
compression (50,000 years), particularly in the neutral mass (Fig.[TO| 
top) and the mean lateral velocity (Fig.[TT|bottom). This is probably 
due to the more vigorous thin-shell instabilities that are seen in the 
higher resolution run. 



3.6 Assessment of the numerical limitations of our 
simulations 

The finite nature of the computational resources available to us 
means that our fixed-grid simulations have a rather modest number 
of grid points in each dimension (127 to 1002). This can potentially 
compromise the reliability of our results both on small scales and 
on large scales. At the smallest scales, we cannot resolve physical 



processes that occur on scales smaller than our cell size. To some 
extent, this is ameliorated by basing the governing equations on 
conservation laws and by careful treatment of the sub-grid physics, 
as in the special case of ionisation fronts ( [Mellema et al.|2006b[ >. 
However, other problems are less easily sidestepped, such as the 
grid viscosity and resistivity, which are many orders of magnitude 
larger than the true viscosity and resistivity of astrophysical plasmas. 
By comparing the behaviour of the same simulation carried out at 
different resolutions, we have determined that the lowest resolution 
that we employ (254 x 127 x 127) is broadly sufficient for modelling 
the global evolution of the globule, even though some of the finer 
scale behaviour is suppressed. In particular, the most compact con- 
figuration of the initial globule implosion is barely resolved in the 
low-resolution simulations and the subsequent fragmentation of the 
neutral globule is delayed and less vigorous. 

At the largest scales, the problem is that our computational 
domain is not many times larger than the size of the globule, so that 
the treatment of the boundaries of the simulation box may seriously 
affect the interior solution. The outflow boundary conditions that we 
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Figure 10. Evolution with time of the total mass of neutral gas (top) and 
ionised gas (bottom). Different line types and symbols show results for 
different model runs, as indicated in the key. 
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Figure 12. As Fig. |10| but showing the evolution with time of the mean 
magnetic field in the neutral gas. Top panel shows x-component of field. 
Bottom panel shows y-component of field. 
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Figure 11. As Fig.[To| but showing the evolution with time of the velocity of 
the centre of mass of the neutral gas. Top panel shows velocity component 
along the .r-axis. Bottom panel shows velocity component along the t/-axis. 
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Figure 13. As Fig. |10| but showing the evolution with time of the mean 
plasma parameter, scaled by the initial value /3o, for neutral gas (top) and 
ionised gas (bottom). Note the different vertical scales in the two graphs. 
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Figure 14. As Fig.|5]but showing the effects of boundaiy conditions on globule evolution. All simulations are similar to S80S but at a lower resolution of 
255 X 127 and each is shown at an evolutionary time of 80,000 years. 



employ (also called non-reflecting or absorbing, e.g., |LeVeque|2002l 
§ 21.8.5) are unproblematic for cases where there is supersonic flow 
outward through the boundary, such as is found in the non-magnetic 
run ZOOL. However, in the strongly magnetised runs, the flow near 
the boundaries is frequently sub-Alfvenic and even subsonic, in 
which case no truly satisfactory treatment of the boundary conditions 
exists. We have attempted to quantify the effects that boundaries 
have on our simulations by carrying out various experiments in 
which the positions of the boundaries or the boundary conditions 
are varied. 

For instance, we have carried out an identical simulation to 
S80L, but with the globule-star system shifted 0.5 pc along the 
positive X-axis. This means that less of the tail is included on the grid 
but that the termination shock of the globule head's photoevaporation 
flow is moved well away from the x = boundary. Some very minor 
differences are seen in the detailed structure of the shocked ribbon 
that forms at late times (§ |3.2.3[ l, but the timing of the recombination 
and the nature of the subsequent fragmentation are identical. 

In another series of experiments, we have carried out 2D simu- 
lations, similar to S80S, but at lower resolution, for three different 
sets of boundary conditions (Fig. |14| l. In the control experiment 
(top panel) the boundary conditions on all faces are the normal 
outflow conditions. In a second experiment (middle panel) we im- 
pose reflection boundary conditions at x = 0, which is equivalent 
to there being two globules along the x axis, one on either side of 
the star at (x, y) = (±0.5, 0.0). In this case, the shock structures in 
the ionised photoevaporation flow from the head of the globule are 
considerably modified from the standard case, since the flow is now 
shocking against the mirror-image flow from the other side of the 
star. The magnetic field lines can no longer be pushed out of the 



computational domain through the x = boundary and the resultant 
increase in magnetic pressure means that the termination shock of 
the photoevaporation flow occurs closer to the globule than in the 
free-outflow case. 

In a third experiment (bottom panel) we have imposed periodic 
boundary conditions aXy = Q and y = y^a.^^, which is equivalent to 
stacking a series of identical globules on top of one another. Again, 
the shock structures in the photoevaporation flow are modified since 
as well as shocking against the ambient medium in the -x direction, 
the flow also shocks against the equivalent flows from the "other" 
globules in the ±y directions. Note that the formation of a dense neu- 
tral shell and shadow zone around (x, y) = (1.0, 0.0) is an unphysical 
artefact of the fact that we have only employed periodic boundary 
conditions for the MHD and not for the radiative transfer. In reality, 
the shell should not form, since the shadow would be ionised by 
the "other" star that is located off" the grid at (x, y) = (0.0, -0.5). In 
this experiment, the accretion flow onto the tail is less dense than 
in the standard case, since there is no longer an infinite reservoir of 
material outside of the box that can be sucked in along field lines 
through the y boundaries. 

Note that despite the changes in the structure of the photoevap- 
oration flow in the experiments with diff'erent boundary conditions, 
the evolution of the neutral globule remains hardly affected. The po- 
tential for "sucking in" material through the boundaries is the least 
satisfactory property of the outflow boundary conditions (which 
could equally well be called inflow conditions). However, we have 
determined that the efi'ects of this on our simulations are gener- 
ally neglible. One exception is the late-time evolution of run S80S 
(Fig.jsj, where, after / = 60, 000 years, a stream of inflowing mate- 
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rial can be seen to develop from the y = l.O boundary at x ^ 1.75, 
which causes distortions in the far tail after t = 100, 000 years. 

The only reliable way of eliminating spurious boundary ef- 
fects would be to make the computational domain so large that all 
physical quantities have become approximately uniform before the 
boundary is reached. However, this approach is computationally 
unfeasible with the uniform-grid algorithms employed in this paper. 
Furthermore, real photoevaporated globules do not exist in isolation, 
but are embedded in complex dynamic environments. Therefore, 
exploring a range of possible boundary conditions, as we have in 
this section, is the best way to determine to what extent the globule 
evolution is affected by its large-scale environment. 



4 COMPARISON WITH OBSERVATIONS 

Magnetic fields have been measured or inferred by various direct 
and indirect methods, both for cometary globules, which show evi- 
dence of interaction with the radiation field of nearby massive stars, 
as well as for dark globules, which show no evidence of such an in- 
teraction. Measurements of both types of globule are relevant to our 
simulations since the dark globules represent the initial conditions 
before the photoevaporation process starts. 

Magnetic field measurements in dark globules are typically 
based on polarisation measurements of the far-infrared emission 
from dust grains (.Henning et al.||2001[ [Vallee et ar]|2003[ |Wolf| 
|et al.|2003[|Vallee & Fiege|2007) , which probes dense gas that is 
very optically thick at visual wavelengths. The projected magnetic 
field direction on the plane of the sky can be determined from 
the orientation of the polarisation vectors, while the field strength 
can be estimated from the dispersion in polarisation direction by 
the Chandrasekhar-Fermi method, so long as the gas density and 
turbulent velocity dispersion are known. The inferred field strengths 
and densities range from 100 to 300 /jG, and from 10^ to 10^ cm"', 
respectively, in globules with total masses from 0.3 to 100 Mq and 
sizes around 0.1 pc. 

The magnetic field in isolated cometary globules has been in- 
vestigated in a series of papers ( Sridharan et al. 1996; Bh att|1999[ 
|Bhattetal.|2004 l via the visual-band polarisation, pv, of background 
stars. As with the far-infrared measurements, this directly reveals the 
projected field direction, and can be also used to estimate the field 
strength from the pv/Ay ratio, where Ay is the visual-band extinc- 
tion, if one assumes that the grain alignment efficiency is the same 
as in the general ISM. In contrast to the far-infrared measurements, 
this technique samples the lower density gas in the globule envelope. 
The only globule where the field strength has been measured by 
this technique is CG 22 in the Gum Nebula jSridharan et al.| 1996b , 
giving a value of 30 fiG in gas with a mean density of 200 cm"'', 
together with some evidence that the field is stronger (70 pG) in the 
higher density (n > 1000 cm"') core. This globule, with a radius 
of 0.6 pc, is considerably larger than the dark globules discussed 
above, but is comparable in mass (~ 15 Mq). The field direction 
is parallel to the globule tail in this case, but in another globule in 
the same region (CG 30-31; pia"tt|1999[ l the field is found to be 
perpendicular to the tail|^ 

In summary, over changes of four orders of magnitude in glob- 
ule density (~ 100 to ~ 10* cm"'), the linear polarization observa- 
tions show only about one order of magnitude variation in globule 



' In a third globule with a measured field direction, CG 12 jfihatt et al.| 
|2004} , the cometary shape is believed to be due interaction with a supernova 
remnant, so it is not directly relevant to the current models. 



magnetic field strengths (6 ~ 40 to ~ 400 pG, after correcting for 
the unobserved line-of-sight component). Assuming a similar gas 
temperature in all the globules, this implies considerable variation 
in the plasma /3 parameter (ratio of thermal pressure to magnetic 
pressure), ranging from y6 ~ 0.01 for the larger, low-density globules 
tOyS ~ 1 for the more compact, dense, dark globules. However, this 
picture is rather different from the results of Zeeman measurements 
of the line-of-sight component of magnetic field in molecular cloud 
cores l |Crutche^| 1999; He iles & Crutcher |2005 Falgaro ne et al.| 
[2008) l, which suggest a roughly constant value of p ~ 0.05 up to 
the highest measured densities (~ 10^ cm"'). It is unclear whether 
this discrepancy is an artefact of the difference in observational 
technique or represents a real difference between the two classes 
of objects. The molecular cloud cores do tend to be much more 
massive (~ 10-1000 Mq) than the isolated globules (~ O.I-IO Mq). 

From our numerical results (§|3|, we find that although devia- 
tions from cylindrical symmetry in the compressed globule occur 
forySo = 0.1, the strongest effects of the magnetic field on the glob- 
ule evolution require Po = 0.01. However, magnetically dominated 
globule evolution may be more common than this fact suggests, 
since classical, isolated cometary globules, such as those in the 
Gum Nebula (Reipurth 1983 1, are typically exposed to a weaker 
ionising radiation field than was employed in our simulations. Two- 
dimensional simulations of globules exposed to a weaker radiation 
field ( |Williams|2007) show that magnetic effects can be important 
even when po ~ I for the case where the ionisation front in the 
ambient medium is D-type when it passes the globule. 

In addition to the clearly isolated globules discussed above, an- 
other class of objects that can be compared with our simulations is a 
whole range of structures seen inside or at the edges of H ii regions, 
which, according to their morphology, may be variously described as 
bright rims I Pottaschj 19 5611, fingers/columns/pillars/elephant trunks 
( |Minkowski| 1 949[ [Osterbrock 1957; Rathborne et al. 2004, here- 
after simply pillars), or bars (O'Dell & Yusef-Zadeh|2000^ . Some of 
the diversity in terminology is spurious, with multiple terms being 
employed for identical or similar structures, sometimes within the 
same paper. Similar objects that have been observed at smaller size 
scales have given rise to yet further profusion of terms: evaporating 
gaseous globules, or EGGS ( [Hester et al.|I996 1, proplyd-like fea- 
tures (De Marco et al.|2006^ , and globulettes l |Gahm et al.|2007[ ). In 
all cases, the underlying physical process is probably the same: a 
dense condensation in the neutral gas slows down the propagation of 
the ionisation front, causing it to wrap around and further compress 
the obstacle. It is usually unclear in individual objects whether the 
dense neutral condensation was pre-existing ( Melle ma et al. 2006a) , 
or whether it formed due to an instability in the propagation of the 
ionisation front ( jGarcia-Segura & Franco|I996[ |Williams|1999[ >. 
Numerical experiments ( Williams et al.|200I^ show that the resul- 
tant structure and evolution of the photoevaporated condensation 
are remarkably similar in the two cases. 

There is one sub-class of these objects that seems to show a real 
physical difference from the others: the linear, bright emission fea- 
tures sometimes referred to as bars, of which the Orion Bright Bar 
is the prototypical example ( |O'Dell|2001[ and references therein). 
Unlike the majority of globules and pillars, these features have a 
long axis that is oriented approximately perpendicular (in projec- 
tion) to the direction of the incident ionising radiation. This is rather 
reminiscent of the emission structures produced in our simulations 
S80H and S80L (Figures|4]and[6](, suggesting that a strong perpen- 
dicular magnetic field may play a role in the formation of bar- like 
features. Magnetic field measurements in the Orion Bar ( [Gonatas] 
let al.|I990[[Houde et al.(2004i|Kusakabe et al.|2008| > are broadly 
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consistent with tiiis scenario: unlike in the main Orion molecular 
filament, where the field orientation is close to the plane of the sky 
and projected along the filament length, the field in the Bar is at 
an angl^to the projected long-axis of the Bar, with indirect evi- 
dence that its orientation is closer to the line of sight. However, the 
true three-dimensional orentations of the Bar, magnetic field, and 
radiation field are far from clear and require further study in order 
to establish whether or not the Orion Bar has been shaped by the 
effects of magnetically modified photoevaporation and implosion. 

The initial conditions employed in the globule simulations 
reported in this paper are highly idealised. In particular, the assump- 
tion of a uniform magnetic field threading the entire simulation 
volume is unlikely to be satisfied in real objects. In a companion 
paper ( [Arthur et al.|2009 |l, we report on global simulations of H ii 
region evolution in highly inhomogeneous, turbulent, magnetised 
media. In those simulations the formation of globule-like photoevap- 
oration flows is an essential feature, as has previously been studied 
in the non-magnetic case ( [Mellema et al.|2006a[|Mac Low et al.| 
[2007l >. 



5 CONCLUSIONS 

We have carried out the first three-dimensional radiation-magneto- 
hydrodynamic simulations of the evolution of a dense, magnetised 
globule that is exposed to a source of ionising radiation. For the 
globule parameters that we employ, we find that when the initial 
magnetic field in the globule is strong (magnetic pressure > 100 
times gas pressure) there are significant effects on the globule evolu- 
tion, such as the following: 

(i) Strong deviations from cylindrical symmetry in the shape 
of the radiatively imploded globule. These deviations are greatest 
when the initial magnetic field orientation is close to perpendicular 
to the direction of the ionising radiation, in which case the shocked 
globule adopts a flattened, plate-like form. In the case of a more 
obliquely oriented initial field, the shocked globule "curls up" to 
form a comma-shaped structure. 

(ii) Modification of the rocket effect, which, in the non-magnetic 
case, accelerates the globule away from the ionising source to reach 
velocities of > 10 km s"'. For a close to perpendicular field orienta- 
tion, the globule acceleration is reversed at late times, whereas for 
an oblique field orientation the globule trajectory is merely deflected 
from a purely radial motion. 

(iii) Magnetic confinement of the ionised photoevaporation flow 
from the globule head, which eventually leads to trapping of the 
ionisation front in a shocked shell, and the temporary recombination 
of parts of the ionised flow. The subsequent evolution is complex, 
but may lead to fragmentation of the shocked shell and the formation 
of secondary mini-globules. 

In the case where the initial magnetic field direction is exactly 
parallel to the direction of the ionising radiation, only the last of 
these three effects applies. In all strong-field cases that we have 
considered, the ionisation front jump conditions on the head of 
the globule are generally of the "extra strong" type, in which the 
magnetic field becomes aligned perpendicular to the ionisation front. 
The densest parts of the ionised photoevaporation flows are always 
found to be gas-pressure dominated, although the ionised gas as a 
whole is magnetically dominated. 

* Note that even if the magnetic field is perpendicular to the Bar's axis, this 
will not generally be true in projection. 



When the magnetic field is weaker (magnetic pressure < 10 
times gas pressure in the initial globule), then the globule evolution is 
more similar to the non-magnetic case, with the principal difference 
being a moderate flattening of the compressed globule. 

We find evidence that magnetic effects may be important in the 
formation of bright, bar-like emission features in H ii regions. 

Our simulations include a careful treatment of the heating and 
cooling in the neutral/molecular gas and we find typical temperatures 
of 50 to 100 K for the densest material, maintained by a combination 
of shocks and x-ray heating. This is warmer than has been assumed 
in previous studies, suggesting that gravitational instability of the 
imploded globule may be less important than has been commonly 
inferred. 
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APPENDIX A: TREATMENT OF GAS HEATING AND 
COOLING IN THE PHAB-C^ CODE 

Al Heating 

The volumetric heating of the gas in our numerical simulations is 
the sum of various terms, each of which represents the effects of a 
dilferent type of radiation: 



H - Hp\]\i + Hy:\]v + Hy + Hm + Hr 



(Al) 



For the gas heating by ionising EUV photons Hevw, we con- 
sider only the photoionisation heating of hydrogen, which is calcu- 
lated in a similar way to the photoionisation rate in equation (jSj of 
§0 



ay{AnJylhv){hv - hvo) dv ergs ' cm ' (A2) 



where /i„ is the number density of neutral hydrogen, cr,, is the pho- 
toionisation cross-section and is the local mean intensity of the 
ionising radiation field, both functions of the photon frequency y, 
as discussed in § |2.1| Note that the hardening of the radiation field 
as one approaches the ionisation front is fully taken into considera- 
tion. In fully ionised gas, //guv dominates the total heating by many 
orders of magnitude. 

In neutral and molecular gas, the ionising radiation is no longer 
present and other heating mechanisms become important. All the 
remaining terms that contribute to the total heating were determined 
by running a large grid of static, constant-density photodissociation- 
region (PDR) models using the plasma physics code Cloudy (Fer- 
|land et al.|1998| l, in which the gas density and incident radiation 
field were varied over a wide range (see Fig. |AI| and § |A3[ >. In all 
these models, the gas-phase abundances and dust properties were 
set at those appropriate to the Orion Nebula. 

For gas close behind the ionisation front, non-ionising far- 
ultraviolet (FUV) radiation with photon energies of around 6 to 
13.6 eV is the most important heating agent, principally via absorp- 
tion by dust grains (including PAHs) and molecular hydrogen lines. 
We find that the resultant heating rate, Hpuv, can be approximated 
as 

1.9x 10-2''«Goe"''"''' 



1 -I- 6A(Go/n)e-^ '>^v 



erg s cm 



(A3) 



where Go is the strength of the FUV radiation field in units of the 
Habing flux (1.2x 10' photons cm"- s"' in the range 912 to 2000 A) 
and Ay is the visual-band optical extinction in magnitudes: Ay = 
l.OSScTvN, in which ctv is the dust extinction c ross-section (assumed 

J nd£ is the 



B aid win et al. 



1991 



and A' : 



column density of hydrogen nucleons. When the gas density, n in 
cm"^, is large compared with the local FUV field, Gqc 



the 



denominator in equation l |A3[ l is close to unity, in which case the 
per-particle heating rate, //puv /" is proportional to the local FUV 
intensity^ For smaller densities than this, the gas-grain thermal 
coupling is inefficient, leading to a reduction in the heating rate. 
FUV heating makes an important contribution in PDRs to a depth of 
Ay = 2-6 depending on the strength of the incident radiation field. 

For a single point source of radiation,the unextinguished FUV 
field is calculated in a similar way as for the ionising radiation in 
eq. l|6j: 

Gfuv 



1.2 X 10' Go 



(A4) 



' Note that a different grain size distribution from the Orion-type grains 
used here would give a steeper decline with Ay. 



where Qfuv is the FUV photon luminosity of the source. For a 
typical O star, Qfuv/Gh is in the range 0.5 to 1.0, depending on 
the spectral type. The column density, N, which is required in or- 
der to determine the local extinction. Ay, is calculated by a short- 
characteristic method in the same way as for the ionising radiation 
( Mellema et al. 2006bj , but using the total hydrogen density in place 
of the neutral hydrogen density. In this treatment, the diffuse non- 
ionising radiation is simply ignored, although it could in principle 



be accounted for by replacing the term Goe in equation (A3 



by a local intensity calculated from a more realistic solution of the 
radiative transfer. 

Hard x rays easily penetrate large columns of gas and can 
therefore be important heating agents deep inside the PDR, in fully 
molecular gas. We have carried out fits to Cloudy models illumi- 
nated by x rays from coUisional ionisation equilibrium plasmas with 
temperatures between 6 x 10* and 4 x 10' K, determining that an 
appropriate approximate heating rate is 



Hx = 6x lO-"nFx 



erg s ' cm \ 



(A5) 



where Fx is the unattenuated x-ray flux, with units erg cm"- s"', in 
the band from 0.5 to 8.0 keV. The Cloudy models extend to Ay = 20 
and show no clear fall-off of Hx with Ay, except for at low gas 
densities. However, we caution against applying equation (|A5j at 
significantly higher columns than this, since attenuation of the x rays 
must become important eventually. 

For sufficiently high gas densities (n> 10^ cm"^) the gas and 
grains are well coupled, so that stellar radiation that is re-processed 
to the infrared at Av ~ 1 and then re-absorbed by dust at larger 
columns will indirectly heat the gas there. From fits to Cloudy 
models, we determine a heating rate for this component of 



HiR = 1.1 X 10-"n(l + nJn)-^Goe-°°^^'' 



erg s cm 



(A6) 



where ni = 3 x 10"* cm"^. 

Finally, when the radiation field is very weak, the baseline 
heating is given by cosmic rays, which we assume to give a constant 
per-particle heating rate such that 



He 



5x10- 



erg s ' cm ^. 



(A7) 



A2 Cooling 

The radiative cooling of the gas is also a sum of various different 
terms: 



(A8) 



In fully ionised gas, collisionally excited optical lines of ionised 
metals dominate the cooling. We approximately model this cooling 
term, Lm+ by multiplying the emission of the typically dominant 
line ([On] 4363 A) by a factor 3 to account for all the other lines 
piroetal.|I995| Appendix A). 



2.905 X 10 



erg s ' cm ^, 



(A9) 



where zo is the oxygen abundance, T[ = 33, 610 K and T2 = 2180 K. 
We ignore any changes in the cooling rate due to second and higher 
ionization of metals (for instance, from to in the interior 
of the ionized region, simply assuming that the strongest cooling 
lines always follow equation ( |A9| l. A similar approach is used for 
collisionally excited lines of neutral metals: 



L^o = 4.477 X I0-^°zon,n„e-'^'''^^''^'''^^' 



erg s ' cm ^, 



(AlO) 



where T3 = 28,390 K and T4 = 1780 K. The product n^n„ is 
proportional to Xi(l-Xi), where x, is the hydrogen ionization fraction. 
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Figure Al. Heating versus Ay for Cloudy models (left panel) and simplified fits (right panel). Solid lines show models including x-ray illumination, dashed 
lines show models with no x rays. Other parameters of the Cloudy models are shown in the key, see text for further details. 
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Figure A2. Cooling versus T from Cloudy models (left panel), together with simplified fits (right panel). 



so that L^o only makes a significant contribution inside ionization 
or recombination fronts, where ~ 0.5. 

We also include cooling terms due to free-free and free-bound 
transitions of ionized hydrogen 

Lh+ = ne«p/H+(7') ergs"'cm"-\ (All) 

and coUisionally excited lines of neutral hydrogen 

Ljjo = ncn„/Ho(r) erg s"' cm"^, (A12) 

where fu*(T) and /noCr) are taken from |Hummer| ( [T994l l. 

For high gas temperatures (above 50,000 K) we adopt the 
following fit to the coUisional ionisation equilibrium cooling curve: 

LciE = 3.485 X 10"'^zor"''"(l - e""""'^'"") erg s"' cm"^ 

(A 13) 

which is mainly due to collisional lines of highly ionized metals (e.g., 
[Dalgamo & McCray|1972^ . Note that this term does not include the 
bremsstrahlung contribution that becomes important at very high 
temperatures, since that is already included in Lh+ . 

For neutral and molecular gas, we adopt a similar approach 
to that described above for the heating, fitting the cooling deter- 
mined from Cloudy PDR models as a function of gas density and 
temperature (Fig. |A2[ l. The resultant fit is 

LpDR = 3.981 X io-2T„^-6T°-Se-ToMIT erg s"' cm"^ (A14) 



where Toin) = 70 + 220(n/10*^)' 



A3 Validation of simplified fits to equilibrium PDR heating 
and cooling 

Figure [Ar| compares the simplified heating functions that we have 
adopted for neutral/molecular gas, //puv + Hx + HlJ^ + Hcr, with 
the results of detailed PDR models calculated with Cloudy. One set 
of Cloudy models was calculated using the following components 
for the incident radiation: 

(i) A 40,000 K black body with Qu = 10"' s"' to represent the 
O star photosphere. 

(ii) An optically thin x-ray source with luminosity Lx,i = 2.29 x 
10^' erg s"' and temperature Txj = 1.6 x 10' K to represent circum- 
stellar emission from the O star i jFeigelson et al.|2005,,Stelzer et al.| 
|2005l l. 

(iii) A second optically thin x-ray source with luminosity Lx,2 = 
3.47x 10^' erg s"' and Tx,2 = 2.5x 10' K to represent chromospheric 
emission from low-mass stars in the accompanying star cluster 
(Flaccomio et al.|2003> . 

(iv) The standard cosmic ray background, equivalent to an H2 
ionization rate of 5 x 10"" s"' ( [Williams et al.|l998t 
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A second set of Cloudy models was calculated with a radiation 
field identical to the above but without any x rays. In both cases, 
the gas phase abundances and grain properties are assumed to be 
those found in the Orion Nebula ( [Baldwin et al.|1991[ |van Hoof] 
[etal. 2004), including a PAH component as in |Abel et 57(2008^. A 
detailed model of molecular hydrogen is included in the calculations, 
as described in ( |Shaw et al.|2005| ). The models are all calculated 
for a plane-parallel geometry, with the total hydrogen density held 
constant at 10^, 10"*, or 10* cm"', with the incident flux calculated 
for distances of r = 0.1, 0.2, 0.5, and 1.0 pc from the ionizing 
star. For components|(i)|and |(ii)[ the flux is calculated assuming a 
point source of radiation, whereas for component |(iii)[ the source is 
assumed to be extended over the stellar cluster or radius = 0.3 pc 
and the flux is approximated as Fx,2 = ix,2/4;r(r^ + r). 

Figure |A1| shows that our approximate heating functions do 
a very good job of reproducing the results of the Cloudy models. 
Note that the lower density Cloudy models include an extended fully 
ionized portion at low Ay (visible as a plateau of high constant heat- 
ing in the graph), in which the heating will be dominated by //euv> 
which is not included in the fits shown in the right panel. Apart from 
in the ionized gas, the region with Ay < 6-10 is dominated by -f/puv^ 
whereas Hx dominates at higher depths for all models that include 
X rays. In the models without x rays, it is Hm that dominates at depth, 
except for at low densities or large distances from the star, where 
HcR becomes important. Additional features are seen in the heating 
curves of the Cloudy models that correspond to the dissociation 
fronts of molecules, particularly H2 and CO. No attempt is made to 
reproduce these low-amplitude features in our fits. Another caveat 
is that the Cloudy models that we have used in constructing our fits 
all assume static equilibrium and do not take into account advective 
and time-dependent effects. Such elfects will be most important at 
discontinuities, where they tend to increase the heating rate with 
respect to the equilibrium calculation ( [Henney et aL|2005] l. In the 
case of the ionization front, this effect is fully taken into account in 
our simulations since the neutral hydrogen density that enters into 
//euv (eq. ( |A2[ )) is calculated from the fully time-dependent equa- 
tion The molecular hydrogen dissociation front, on the other 
hand, is not explicitly included in our fits and advective effects there 
are not accounted for. However, for the case of photodissociation 
regions illuminated by O stars, the changes introduced by advection 
are expected to be rather small (Stoerzer & H ollenbach| 1 998 1 . A 
notable exception to this will occur in cases where the ionization 
front and dissociation front merge, which tends to occur for low 
ionization parameters and hard incident spectra ( jBertoldi & Draine] 
[1996). In such cases, a much more careful treatment of the heating 
is necessary, as in |Henney et al.|p007[ l. 

Figure |A2]presents a similar comparison for our approximate 
cooling function Lpdr. It can be seen that this function does a good 
job of reproducing the principal variation of the cooling with den- 
sity and temperature. At temperatures below ~ 100 K, the cooling 
is prinicipally due to millimetre CO lines (at higher densities) or 
far-infrared neutral C lines (at lower densities), while at warmer tem- 
peratures it is dominated by far-infrared neutral O lines (at higher 
densities) or singly ionised C lines (at lower densities). At temper- 
atures greater that 3000 K, the cooling begins to be dominated by 
the optical and UV lines that we include in Lyo and L^^o using the 
non-equilibrium time-dependent ionization fractions, so this portion 
of the cooling is not included in the equilibrium fits shown in the 
right panel. At the lowest temperatures, one sees that the cooling 
in the Cloudy models is no longer uniquely determined by T and 
n, especially for the lower densities. Most of this extra variation, 
which we make no attempt to reproduce in our fits, is due to changes 



in the molecular abundances. Also shown in the figure is the curve 
of Koyama & Inutsukal ( |2002[ >, as corrected by |Vazquez-Semadeni| 
|et al. 1 2007[ l, which has been frequently employed in numerical 
simulations: 



Lain, T) = 2x lO^'V [e 



-1. 184x10' /(r+IOOO) 



+ 1.4 X lO^'T'/^e"'-'^) erg s"' cm"' (A15) 

It can be seen that this curve is a reasonable approximation to 
the cooling at low densities, but that it seriously overestimates the 
cooling rate for n > 10' cm"', being 100 times too large for T = 
lOOK, «= 10* cm-'. 

In summary, we believe that the heating and cooling functions 
that we have introduced in this appendix represent a substantial 
improvement over what has been commonly used in the previous 
literature, which has generally either assumed an isothermal equa- 
tion of state for neutral/molecular gas (e.g.,|Esquivel & Raga 2007]l 
or used H = Hqr and L = Lki (e.g., |Krumholz et al.|2007p . The 
typical neglect of FUV/optical dust heating (important for columns 
with Av < 5), of x-ray heating (important within ^ 3 pc of a typical 
O star) and of collisional deexcitation of the principal cooling lines 
(important for « > 10"* cm"') will tend to result in an underestimate 
of the temperature of shocked neutral/molecular gas. Since we have 
broken down the heating function into terms due to different wave- 
length bands, our results can be used in the case of illumination 
by stars of different effective temperatures and with different x-ray 
luminosities. 
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